Nucleation of cracks in a brittle sheet 
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We use molecular dynamics to study the nucleation of cracks in a two dimensional material 
without pre-existing cracks. We study models with zero and non-zero shear modulus. In both 
situations the time required for crack formation obeys an Arrhenius law, from which the energy 
barrier and pre-factor are extracted for different system sizes. For large systems, the characteristic 
time of rupture is found to decrease with system size, in agreement with classical Weibull theory. 
In the case of zero shear modulus, the energy opposing rupture is identified with the breakage of a 
single atomic layer. In the case of non-zero shear modulus, thermally activated fracture can only 
be studied within a reasonable time at very high strains. In this case the energy barrier involves 
the stretching of bonds within several layers, accounting for a much higher barrier compared to the 
zero shear modulus case. This barrier is understood within adiabatic simulations. 



I. INTRODUCTION 

While our current understanding of fracture begins 
with the ideas of Griffith in 1921 [1 , the study of its 
atomic mechanism has attracted a large amount of at- 
tention in recent years. For example, corrections to Grif- 
fith's results for a crack in a brittle material have been 
proposed and verified with atomistic simulations [3J |31 13] • 
Also, large scale simulations have been used to study 
dynamical fracture [5] El IZj- This increase in interest 
in fracture is partly due to computer simulations which 
promise an understanding at the atomic level description 
of the phenomena. However simulations face a funda- 
mental problem [51 [5]: many atomic deformations are 
thermally activated and therefore involve long timescales 
which are difficult to simulate. 

Most simulations overcome this problem by studying 
fracture with a pre-existing crack. In that case, crack 
growth is a driven phenomena and there is no energy bar- 
rier to be overcome. Only a few simulations have been 
used to study the formation of cracks at non-zero finite 
temperature without pre-existing cracks: void formation 
has been observed in 3D simulations of strained binary 
Lennard- Jones systems |10j . and simulations for the rate 
of crack nucleation have been performed in a 2D spring 
network Experimentally, the rate of crack nucle- 

ation in heterogeneous materials has been found to obey 
an Arrhenius law with an energy barrier scaling accord- 
ing to Griffith's results [HHS]. 

In the present work, we address the nucleation of cracks 
in a brittle two-dimensional material, i.e. a sheet with 
a thickness of one atomic layer, through Langevin dy- 
namics. The rate constant for the nucleation of cracks 
follows an Arrhenius law, from which the energy bar- 
rier is extracted. Two situations in a square lattice are 
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studied: atomic interactions restricted to first neighbours 
and interactions extended to second neighbours. In the 
former case, the shear modulus of the solid is zero and 
the energy barrier is shown to be independent of system 
size: the breakage of a single bond propagates to the rest 
of the solid without any cost. In the latter case the shear 
modulus is non-zero and a finite size crack has to nucle- 
ate before rupture can propagate throughout the system. 
These two situations will be referred to as chain-like and 
solid-like models, respectively. 

This paper is organized as follows. In the next section 
we describe the models used in this paper followed by a 
description of how simulations are carried out. In section 
IV we present the results for the chain-like and solid-like 
models. The latter is physically more relevant to describe 
brittle materials and we discuss its energy barrier in the 
context of Griffith theory in section [V] Our conclusions 
finalize the paper. 



II. MODEL 

A stretched one dimensional chain has been previ- 
ously used as a simple model for breakage of poly- 
mers [131113 HnHni nil [H]. Here we extend this model to 
study fracture in 2D brittle solids by bonding the chains 
to each other such as to form a square lattice - see Fig. [T] 
We study samples containing M chains which are made 
of = 100 atoms each. Those chains are stretched in 
the horizontal direction: their constant length is N{a+s), 
where a is the equilibrium bond length and s is the ap- 
plied strain. By constraining the dynamics of atoms 
along the applied strain the system cannot form topolog- 
ical defects and the only mechanism for stress relaxation 
is fracture. Also, by choosing the constraint along the 
applied strain, we expect to be sampling the meaningful 
pathway for fracture while speeding up the simulation 
time. This is verified later in section |IV| where we per- 
form one set of simulations without this constraint. 
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FIG. 1: (Color Online) Schematic representation of system 
containing M x N = 5 x 8 atoms. Atoms are only allowed to 
move in the horizontal direction. First-neighbor interactions 
are represented by straight lines while dotted lines represent 
second-neighbor interactions. Periodic boundary conditions 
are represented by dashed lines. 

A square lattice can be made isotropic by choosing 
the elastic constant between first neighbors to be twice 
as large as the elastic constant between second neigh- 
bors f2(r. To fulfil this condition we chose the following 
form for the Lennard- Jones potential between first neigh- 
bors: 

Vf{r)^e[{a/rY'~2{a/r)% (1) 

and 

K(r)=4e,[(a,/r)i2-(a,/r)«], (2) 

for second neighbors. We use e — 1 and = 36/228e 
for the binding energies and a = 1 and = ^/2a for the 
equilibrium lengths. 

The dynamics of this system are obtained by solving 
a set of Langevin equations for the position Xij of each 
atom: 

'^'^^^= T,kdPixi,3^Xk,i)~V^ + kj{t) (3) 

where F(x) is the force computed from the potential, m 
is the atomic mass and rj is the friction coefficient. The 
random force fi,j{t) is related to rj by the fluctuation- 
dissipation theorem. 

Periodic boundary conditions in the horizontal direc- 
tion imply xoj — xnj and x^+i.j — xij for all j. 
Periodicity is also imposed in the vertical direction to 
ensure that all chains are equivalent: Xi^ = x^^m and 
Xi^M+i = Xi^i. For simplicity we use reduced units. En- 
ergy is given in units of e, distance is given in terms of a, 
and time is given in units of the smallest phonon oscil- 
lation period P = 2it / {12\/{2e/raa?)) of an intact chain 
[14] . Mass is written in terms of m and the friction coef- 
ficient is tuned to ry = 0.25(27r/P). 

Initially all the horizontal bonds have the same length 
a+s and all vertical bonds are at their equilibrium length 
a. The velocity of each atom is chosen randomly accord- 
ing the Boltzmann distribution. The dynamics of the 
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FIG. 2: (UPPER PANEL) Dependence of the potential en- 
ergy for the solid-like model defined by M x A'' = 60 x 100, 
S — 0.065 and T — 0.016. Arrows indicate the instances at 
which the atomic configurations in panels A,B,C and D are 
shown. 



system are obtained by solving numerically the set of 
equations |3] using the velocity- verlet algorithm until 
the solid ruptures. 



III. SIMULATION 

In Fig. |2] (upper panel) we show the time dependence 
of the potential energy for the solid-like model. This en- 
ergy fluctuates around its initial value showing that the 
initial stretched state is a configuration at a local en- 
ergy minimum. Rupture occurs at about 2750 units of 
time when the energy of the system drops abruptly. This 
shows that the system is been driven towards an equilib- 
rium state with lower energy. In panels A,B,C and D of 
this figure, we shown atomic configurations at different 
instances along fracture. In those panels, the incipient 
crack is seen to propagate perpendicularly to the direc- 
tion of applied strain. This indicates that we can use 
the sum of bond lengths along the pathway where frac- 
ture is taking place as an order parameter (j) for fracture. 
For convenience, the sum is taken over all the largest 
bonds percolating vertically along the sample and only 
the horizontal bonds are considered in the sum. Thus, 
initially = M{a + s) and increases until two surfaces 
are formed. 

Notice that the potential energy in the upper panel of 
Fig.[2]shows no apparent precursor behavior for rupture. 
Also, the energy barrier that the system has to overcome 
for rupture to proceed is smaller than fluctuations in the 
total potential energy. Thus it cannot be easily extracted 
from an analysis of the potential energy. To obtain this 
barrier we study the kinetics of the system as it proceeds 
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FIG. 3: Dependence of the characteristic time on the cut-ofT 
value of the order parameter for the sohd-like model defined 
by M X iV = 60 X 100, S = 0.065 and T = 0.016. The 
dashed line separates reversible (left side) from irreversible 
(right side) rupture. 



towards rupture. In particular we measure the charac- 
teristic time of rupture and analyze this quantity from 
the point of view of thermally activated systems. 

To compute the characteristic time of rupture r we 
use an ensemble containing Sq = 1000 samples for the 
chain-like model and 5*0 = 500 samples for the solid-like 
model. Those samples differ from each other by the sets 
of initial velocities and random forces fi.j{t). We will 
need to choose a value of the order parameter cp which 
we will associate with irreversible rupture. The charac- 
teristic time for the incipient crack to reach a particular 
size is computed by tracking the number of chains S^cj), t) 
whose order parameter has not yet reached the value (f> 
at time t. For a fixed value of (j) this function decreases 
exponentially with time, S{(j),t) = Soexp(—t/T{4>)). The 
characteristic time t(0) depends on cj) and is obtained 
from a fit of S{4>, t) to the numerical data. 

The time that characterizes rupture in an irreversible 
manner depends on the arbitrary choice of the cutoff 
value (pc which we associate with rupture. To chose this 
cut-off we show in Fig. |3] the dependence of r on (j>. Two 
distinct regimes are apparent. The first regime occurs 
when (j) is smaller than ~ 65.8 (in units of a). In this 
regime, 4> increases very slowly with time. The under- 
lying physics of this regime is the competition between 
thermal fluctuations, which are responsible for increasing 
crack length, and the restoring force on the atomic bonds. 
The second regime occurs when the order parameter (j) is 
greater than 65.8. Here, t(0) has reached a plateau, and 
(j) increases very rapidly with time. Stress relief of the 
material's bulk is the driving force of this regime which 
requires a larger crack and therefore produces a fast in- 
crease in (p: irreversible rupture has occurred. So, from 



FIG. 4: (Color online) CHAIN-LIKE MODEL - (a) Depen- 
dence of Inij) on the inverse of temperature for different sys- 
tem sizes M . (b) Dependence of tq on system size, (c) De- 
pendence of the energy barrier Ef, on system size. 

Fig.[3]we can determine the value of for which rupture 
becomes irreversible. This value is (pc = 65.8. 

IV. RESULTS 

The nucleation of cracks can be thermally activated 
O [9l [11] such that their occurrence is typical of an Ar- 
rhenius process. Mathematically the characteristic time 
of rupture r, the inverse of the nucleation rate, reads: 

r = roexp(£;fa/fcbr) (4) 

where k^T is thermal energy, Eb is the energy barrier 
the system has to overcome, and tq is the inverse of the 
attempt frequency to rupture. The attempt frequency 
depends on the vibration frequency of the system in the 
metastable wells of the energy landscape through the lo- 
cal curvature of the energy surface [23J [H]. It also de- 
pends on the friction coefficient in the Langevin equa- 
tion [22] . In this section we study the dependence of t 
on temperature for the chain-like and solid-like models 
to extract both Ef, and tq, which are intrinsic quantities 
of the system being studied. 

A. Chain-like model 

In Fig. [4] we show the temperature dependence of t 
for different system sizes M. The strain of this system 
is set to s = 0.05. For each system size, r increases ex- 
ponentially with 1/kbT - in agreement with Eqn.[4j The 
energy barrier E}, and the pre-factor tq were extracted 
from fits to those results. Changing system size strongly 
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affects the pre-factor but has no effect, within error bars, 
on the energy barrier - as can be seen in Figures 
and (c) respectively. 

The dependence of the pre-factor on system-size can be 
understood qualitatively within the scope of a nucleation 
theory for fracture. Intuitively, tq is proportional to the 
average time between two consecutive attempts to nucle- 
ate a crack. Assuming that cracks can nucleate in each 
of the M chains of the system, then tq = l/M [26ll27] . 
Fig. shows the good agreement of simulation with 
this inverse relation. 

The energy barrier can be understood quantitatively 
by assuming that parallel chains are independent from 
each other. Under this assumption, the energetic cost E 
of elongating one atomic bond in a single layer is given 
by [Ml [ig Hi]: 
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E{^)^V{a + s + cj,) + {N-l)v(^a + s-j^^^^) (5) 

where (j) is the deviation of the broken bond length from 
its strained elongation and V{x) is the potential energy 
of an atomic bond. Equation [s] corresponds to the sum 
of potential energy of all the bonds in the layer precur- 
sor of fracture. This equation considers that while one 
of the bonds increases towards rupture by an amount 0, 
the other bonds of the same layer relax by an amount 
4>/{N — 1) towards their equilibrium value. For the pa- 
rameters used in Fig. |4] i.e. TV — 100 and s — 0.05, 
Equation [5] predicts an energy barrier of 0.0564. This is 
in good agreement with our simulations where the barrier 
is approximately 0.054 for all system sizes. 

Independence of parallel chains is the key assumption 
to explain fracture in the chain-like model. This assump- 
tion can be understood as follow. Since the shear mod- 
ulus of this model is zero, no energetic cost is associ- 
ated with sheared configurations in the linear regime. 
Therefore, an individual atomic layer can proceed to- 
wards fracture independently of neighbouring layers un- 
til non- linear effects become relevant. The energy barrier 
opposing this process is related to the cost of increasing 
the length of one of the bonds of the layer - indepen- 
dently of other layers. Only when this bond becomes 
large enough, neighbouring layers are driven towards rup- 
ture in a domino- like process. 

Eqn. |5] results from the competition between the en- 
ergetic cost of extending one bond length of the chain 
and the energetic gain of relaxing the remaining bonds. 
This contrasts with Griffith's calculation where the bar- 
rier is related to the cost of creating more surface and the 
energetic gain due to reducing the strain in the bulk of 
the material. Thus, despite its use in the literature, the 
square lattice with only first neighbour interactions is a 
poor model for fracture in a solid and Griffith's theory 
does not apply to this system. 



FIG. 5: (Color online) SOLID-LIKE MODEL, (a) Depen- 
dence of /n(r) on the inverse of temperature for different 
system sizes M. (b) Dependence of tq on system size, (c) 
Dependence of the energy barrier Ei, on system size. The re- 
sults for the "full 2D" system correspond to a dynamically 
unconstrained model - see text for more details. 



B. Solid-like model 

We now study the solid-like model. In Fig. [5] we show 
the temperature dependence of r for different system 
sizes M and applied strain s = 0.065. As in the pre- 
vious model we study this system by fitting the time of 
fracture to Eqn. [4] obtaining the energy barrier Ei, and 
the pre-factor tq for each system size. Those results are 
shown in Fig. |5]^b-c). 

The pre-factor. Fig. ^h), presents two regimes: for 
systems containing less than 15 chains, i.e. M < 15, tq 
increases with system size; however for Af > 15, the pre- 
factor decreases as system size increases. Those behav- 
iors are related to finite size effects. When M < 15, the 
relaxation region around the crack is of the same size as 
the system. On the other hand, increasing the size of the 
solid above 15 layers implies that more nucleation sites 
are available for rupture, and tq decreases with M, as 
in the previous model [521 HZ). A quantitative assess- 
ment of the pre-factor would involve the generalisation 
of Kramer's result to higher dimensions [28 . This cal- 
culation was performed successfully to study rupture in 
a one-dimensional chain [14J but its application to the 
present model is beyond the scope of this paper. 

In Fig. ^c) we show the energetic cost for nucleating 
a crack in the solid-like model as a function of system 
size M. For solids smaller than M — 15, the energy bar- 
rier increases considerably with system size: more than 
150% in Fig. ^c); while for solids larger than M = 15, 
the increase is only marginal and show a saturation trend 
at El, ~ 0.14 - indicating that finite-size effects become 
negligible. This value is comparable to the adiabatic bar- 
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rier. This barrier is computed by extending the length 
of one bond in small steps and restraining it while the 
other bonds are relaxed at zero temperature. In this 
process the energy increases until the critical crack is 
formed. The maximum energy seen in this process cor- 
responds to the energy required to nucleate the crack 
at zero temperature. This adiabatic energy for the dif- 
ferent system sizes are represented by squares in Fig. [5] 
Notice that the barrier obtained in our simulations is 
smaller than the adiabatic energy barrier by about 23 %. 
A smaller simulated barrier compared to the adiabatic 
case has also been observed for one dimensional systems 
[IB] . A possible explanation for this discrepancy might 
be that a zero temperature calculation does not account 
for entropy which plays a role in the free energy opposing 
rupture in system with multiple degrees of freedom. 

One important simplification imposed in our model 
with respect to two-dimensional solids consists in con- 
straining the dynamics of atoms to one dimension. How- 
ever by imposing this constraint along the direction of 
applied stress, we expect to be sampling the meaningful 
pathway for rupture of a 2D-solid. To verify this state- 
ment we performed a set of simulations on a M x A'^ = 
50 X 100 system where the constraint on the motion of 
atoms was removed. The results of those simulations are 
shown in Fig. [5][b-c) and are referred to as full 2D. No- 
tice that the full 2D system has a much lower pre-factor 
than our constrained system. A discrepancy in the pre- 
factor is expected since it is related to the vibration of 
the system, and therefore its dynamics, which is differ- 
ent in both models. However the energy barrier of our 
constraint model and the full 2D are equal within error 
bars. We are therefore confident that our constrained 
model can be used to study the energetic behavior of 2D 
solids. In the next section we discuss the simulated en- 
ergy barrier in the context of Griffith theory for rupture 
and adiabatic simulations. 



V. SOLID-LIKE MODEL AND GRIFFITH 

The introduction of a crack of size L in a solid charac- 
terised by a Young's modulus E and subjected to a stress 
a will result in a stress-energy relief of na^L^ /2E. But 
this crack will also involve a cost of 27L, where 7 is a 
surface energy, such that the dependence of the energy 
on the crack size L is [13]: 

EGiL)^^^^ + 2^L. (6) 

This potential energy reaches a maximum when 
dEG{L)/dL — 0. This occurs at the critical value 
Lg = 2i?7/(7rcr^). Beyond this crack length, the crack 
propagates spontaneously to reduce the bulk strain in 
the material until the solid is broken in two pieces. 
The barrier for crack nucleation occurs at this critical 
length: Eg{Lg) = ^ or The Young modulus of 

the solid- like model is E = 77.91 (in units of /a?) and 
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FIG. 6: Relative difference in energy barrier as computed 
from Griffith calculation {Eg) and the adiabatic simulation 
{Ea): {Eg - Ea,)/EG. INSET - dependence of the adiabatic 
barrier on strain in a log-log scale. The linear fit of systems 
with a strain smaller than 0.04 is shown. 



the energy necessary to create two surfaces 27 is equal to 
the energy of two weak and one strong bond per inter- 
atomic distance: 7 = 0.6587. For a strain s = 0.065, the 
Griffith energy barrier is Eg = 0.837, that is, approxi- 
mately six times the value obtained from our statistical 
simulation. 

It is of no surprise that Griffith's calculation is not valid 
for large strains. First, because in this highly stretched 
regime, linear elasticity theory is not valid. Second, for 
strains larger than 0.04 only one atomic layer needs to 
be completely broken in order to initiate the rupture pro- 
cess. In other words, while in Griffith's regime the mech- 
anism behind rupture is the competition between the for- 
mation of new surface and stress relaxation, the physics 
of rupture in the highly stressed regime is the competi- 
tion between stress relaxation and bond stretching at the 
formation of the incipient crack. 

To understand the range of validity of Griffith calcu- 
lation, we performed adiabatic relaxation (T = 0) where 
bonds were mathematically cut within a line perpendic- 
ular to the applied strain. After cutting those bonds, 
atoms were relaxed until the force on each one of them 
was smaller than 1 x 10^^. To avoid finite-size effects we 
increased the size of the system from N=100 at S=0.04 
to N=600 at S=0.01. The relative energy barrier com- 
puted from this process with respect to Griffith barrier 
is shown in Fig. [6] This Figure shows that for strains 
below 0.04, the adiabatic barrier differs at most by 30 % 
from the Griffith barrier. On the other hand, for strains 
beyond 0.04 those barriers are several times different and 
this difference increases with strain. In the inset of Fig. 
|6] we also show in a log-log scale the dependence of the 
energy barrier on strain. For strains below 0.04, those 
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quantities scale with an exponent of -2.26. This is very 
close to the exponent predicted by Griffith's theory: -2. 
The behavior at higher strains deviate from this scaling. 
This clearly shows that for our model Griffith's theory 
is valid for strains smaller than 0.04. Wc note that the 
calculation of the energy barrier with strain dependent 
Young modulus and surface tension - as described in Ref . 
[3] - gave results in greater disagreement than the ones 
of Griffith. 



VI. CONCLUSION 

Due to its simple dynamics and large system size, our 
atomistic simulation of the nucleation of cracks in thin 
brittle sheets is an ideal system for the study of noise 
activated processes and nucleation theory. In particular, 
we found that the energy barrier for crack nucleation in 



a square lattice with only first-neighbor interactions is 
comparable to the barrier of one-dimensional chains due 
to the zero shear modulus of this system. For the more 
interesting case where second-neighbor interactions are 
incorporated into the model, we found an agreement be- 
tween the simulated energy barrier at high strains and the 
one computed from an adiabatic relaxation. This barrier 
involves several layers, accounting for a much higher bar- 
rier compared to the case of isolated chains. We believe 
that extensions of the present study such as to investi- 
gate nucleation of pre-existing cracks would be a valuable 
contribution to understand fracture. 
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